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Abstract Marine vertebrate strandings offer an opportunistic sampling scheme that can 
provide abundant data over long periods. Because the stranding process involves biolog- 
ical, physical and sociological parameters, confounding complicates the interpretation of 
results. The statistical analysis of these data relies on generalized linear or additive models 
in order to infer long-term trends, but does not easily account for drift or variation in 
reporting rates. Here, we capitalized on county-level (administrative) variation following 
the passing of a law for compulsory reporting of stranded marine mammals in France to 
investigate how variation in reporting rates may affect the observed trend in stranded small 
delphinids in the Bay of Biscay. Using a time-series spanning more than 30 years across 
eight administrative counties, we built variance partitioning models for the analysis of 
count data. We discussed the choice of an appropriate likelihood to conclude the Negative 
Binomial useful and interpretable in the context of small delphinid strandings. We 
expanded the model with a recent methodology to detect structural breaks in the time 
series, focusing on overdispersion. We performed statistical robustness checks with respect 
to variations in reporting rates and discuss their causal interpretation in the context of 
observational data. Stranding frequencies increased on average 7-fold over 30 years. We 
conclude that reporting rates to the French stranding network have been stable since the 
early 1990s, and the average 3-fold increase in stranded small delphinids observed in the 
Bay of Biscay since 1990 is due to other factors, including bycatch. Codes and data are 
available to replicate the analysis to other national stranding networks. 
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Introduction 

The European Union issued a Marine Strategy Framework Directive (MSFD hereafter) to 
protect and restore the ecological integrity of marine ecosystems within member states. 
The ambitious goal of the MSFD is to achieve "Good Environmental Status" (GES 
hereafter) of marine ecosystems on a suite of 1 1 descriptors, including biodiversity. GES is 
defined as "the environmental status of marine waters where these provide ecologically 
diverse and dynamic oceans and seas [. . .], and the use of the marine environment is at a 
level that is sustainable [...]" (Article 3.5 of Directive 2008/56/EC of the European 
Parliament and of the Council of 17 June 2008). Among the threats looming on marine 
biodiversity, bycatch or the capture or entanglement of non-target species in fishing gear, is 
foremost (DeMaster et al. 2001; Reeves et al. 2013). Monitoring marine species is para- 
mount to diagnose impediments or progress towards GES. Monitoring marine mammals is 
however challenging (Evans and Hammond 2004). While dedicated surveys are needed to 
obtain crucial parameters such as absolute abundances (Hammond et al. 2013), the low 
temporal frequency at which such surveys can be carried out results in low statistical power 
to detect declines (Taylor et al. 2007). Yet the MSFD aims at monitoring relative progress 
and diagnosing problems at a finer temporal resolution: complementary data are needed in 
order to implement conservation measures in a timely fashion. 

Marine vertebrate strandings offer an opportunistic sampling scheme that is valuable to 
study abundant, rare and cryptic species. Dedicated at-sea surveys for such species are 
logistically difficult and costly (Evans and Hammond 2004; Thompson et al. 2012). 
Stranded cetaceans provide raw material for many scientific purposes: for example pop- 
ulation genetics (Amaral et al. 2012; Bilgmann et al. 2011), reconstructing individual 
reproductive history (Dabin et al. 2008), foraging ecology (Spitz et al. 2012; Dunshea et al. 
2013), documenting temporal trends (Fruet et al. 2012; Pikesley et al. 2012; Truchon et al. 
2013; Williams et al. 2011); explaining extreme stranding events (Fernandez et al. 2012; 
Jepson et al. 2013; Rubio-Guerri et al. 2013; Williams et al. 2011; Wright et al. 2013); or 
estimating at-sea mortality (de Boer et al. 2012; Koch et al. 2013; Peltier et al. 2013; 
Williams et al. 2011). 

Stranding is the result of at-sea mortality, buoyancy, local drift and prevailing winds, 
shore substrate and detection probability (Peltier et al. 2012, 2013; Williams et al. 2011). 
Peltier et al. (2014) described stranding data as 3-dimensional: (1) a biological dimension 
which includes mortality rate and abundance; (2) a physical dimension linked to the drift 
process; and (3) a societal dimension encompassing detection and reporting of stranded 
carcasses by citizens. Interpreting stranding records is difficult as the outcome of interest, a 
standardized number of stranded carcasses, results from many processes interacting along 
these three dimensions (Peltier et al. 2012; Williams et al. 2011). In particular stranding 
probabilities are difficult to estimate in situ, but substantial progress has been made with 
drifter experiments (Koch et al. 2013) or physical modelling of drift and winds (Peltier 
et al. 2012, 2013, 2014). 

In monitoring stranded cetaceans, interest also lies in elucidating the causes behind 
observed patterns. Causal inference usually follows from evidencing a statistically sig- 
nificant relationship over time. Implicit here is the view of causation as "robust depen- 
dence" (Goldthrope 2001). This view amounts to little more than statistics because no 
explicit reference to theory needs to be made (Goldthrope 2001). While monitoring may be 
of political value, it can be scientifically inefficient without the ability to manipulate 
(Yoccoz et al. 2001). The view that only manipulation legitimates causal inference falls 
under the view of "causation as consequential manipulation" (Goldthrope 2001). A causal 
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model goes beyond mere description: it provides means to predict how a response variable 
will change when one or more causal variables are manipulated (Rubin 2006). With 
stranding records, 'time' is not a variable that can be manipulated, nor is it usually a well- 
defined treatment (Cox 1992; Holker et al. 2007). However, a limitation of this view of 
"causation as consequential manipulation" is its lack of reference to a generative process 
operating at a level below that of the observed pattern (Cox 1992; Goldthrope 2001). This 
realization invites for another view of causation: "causation as generative process" 
(Goldthrope 2001). Under this view, both theory (at a proximate and mechanistic level) 
and predictability (at a larger scale) are considered in tandem with robust statistical 
dependence to perform causal inference with observational data. 

Establishing causation requires an interpretable model with good predictive ability, and 
that potential confounders are taken into account to obtain robust results. Because the 
number of stranded carcasses is the result of many underlying processes (Peltier et al. 
2012), any detected trend in stranding records may be due, for example, to varying dis- 
covery and reporting rates (hereafter detection probability) of carcasses over time. Albeit 
direct manipulation is usually impossible, "quasi-experiments" may sometimes seren- 
dipitously arise in observational studies: any event unrelated to the outcome of interest that 
generates variation provides an opportunity to mimic an experimental situation. In par- 
ticular geographical variation may provide opportunities to attempt causal inference with 
observational data. In large countries, the subnational organization of administrative ter- 
ritories offer means to investigate patterns of variations while still holding important 
factors, such as legislation, constant (Barr et al. 2012). We propose to analyze French 
stranding data at the administrative county level to capitalize on known differences in local 
conditions between counties. Eight French counties border the Bay of Biscay, with "Les 
Landes" county differing in one major respect: since 1991, monitoring has been institu- 
tionalized. In this county, local permanent staff is in charge of monitoring and cleaning 
beaches, for aesthetic reasons linked to tourism. We expect reporting rates of stranded 
cetaceans to have remained fairly constant since the 1990s there. 

We built within a Bayesian framework statistical models for analyzing data collected 
over > 30 years of daily monitoring of stranded carcasses of Common dolphins (Delphis 
delphis) and Striped dolphins {Stenella coeruleoalba) along more than 3,000 km of 
shorelines in the Bay of Biscay (Fig. 1). Emphasis was first put on model predictive ability: 
we used statistical models and assess their ability to generate data similar to the observed 
data ("causation as generative process"). Secondly, we checked the robustness of our 
results to the inclusion of variables related to reporting rates. Throughout, we focused on 
modelling overdispersion. Finally, we discussed the biological interpretation of this 
parameter in the context of small delphinid strandings. 

Materials and methods 

Ethics statement 

The study is entirely based on data collected from cetacean carcasses found stranded along 
the French coasts and did not involve observation or experimentation on captive animals 
by any means, nor did it rely on field observation of live animals. 

Observatoire PELAGIS (formerly "Centre de Recherche sur les Mammiferes Marins") 
is the institution permanently in charge of coordinating the French marine mammal 
stranding network ("Reseau National d'Echouage", hereafter RNE) under the decree of 
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Fig. 1 Bay of Biscay. The eight French administrative counties ("Departements") bordering the Atlantic 
Ocean and included in the analysis are emphasized with color. On a South-North gradient, these counties 
are "Pyrenees Atlantiques", "Les Landes", "Gironde", "Charente Maritime", "Vendee", "Loire 
Atlantique", "Morbihan" and "Finistere". The French Exclusive Economic Zone extending in the Atlantic 
Ocean is delineated with a black line 

November the tenth of 2010, jointly taken by the Ministery in charge of the Environment 
and the Ministery in charge of Fisheries, regarding the use of biological data and samples 
collected on stranded marine mammals for scientific research and monitoring purposes. All 
species of marine mammals occurring in waters under French juridiction are protected by 
the decree of July the first of 201 1. 1 Since the early 2000s, the French Ministry of the 
Environment, Sustainable Development and Ecology mandates Observatoire PELAGIS 
(UMS 3462 du CNRS, Universite de La Rochelle) to centralize information, biological 
samples and data on marine mammals in mainland France. It authorizes Observatoire 
PELAGIS staff and trained volunteers (that is, holders of a "Green Card" delivered upon 
completion of a training program dispensed by Observatoire PELAGIS) to manipulate and 
sample cetacean carcasses for scientific purposes. In the present study, sampling only 
involves species determination. This work was carried out following European regulations 
on the use of stranded dead cetaceans for scientific purposes. 

Sample collection 

The RNE is dedicated to the monitoring of marine mammal populations. Around 260 
trained volunteers are currently taking an active part in the network. These volunteers 
make the complete coverage of French coastlines possible. Standardized training of vol- 
unteers by permanent Observatoire PELAGIS staff, which takes place twice a year, ensure 
the homogeneity, comparability and standardization of data collection procedures in the 



1 http://www.legifrance.gouv.fr/affichTexte.do ?cidTexte=JORFTEXT000024396902. 
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field. Data are further checked for consistency or errors before their incorporation in a 
national database managed by Observatoire PELAGIS. Each stranded animal corresponds 
to a datum in the database. 

The RNE was established in the early 1970s, but its modus operandi is considered 
unchanged since the early 1980s. The mandatory reporting of cetacean strandings to RNE by 
town and village administrations (municipalities) was enforced by law in 1988. Among the 
> 36,000 French municipalities, 800 have a geographic access to the sea, and may report a 
stranded animal. We used data collected in the eight counties ("departements") bordering the 
Bay of Biscay (Fig. 1). These counties encompass approximatively 300 coastal municipal- 
ities, all of which have reported at least one stranded cetacean to RNE since its inception. 

Local citizens dwelling in a coastal municipality may spot a stranded cetacean and 
either contact their local authorities or Observatoire PELAGIS directly, which will sub- 
sequently send a RNE volunteer to collect data on the stranded cetacean. Thus densely 
inhabited coastal towns may have higher detection rate of strandings. Demographic data on 
these municipalities have been collected by "Institut National de la Statistique et des 
Etudes Economiques" approximatively every ten year since the 1970s. We downloaded 
demographic data from http://www.insee.fr in order to incorporate this covariate as a 
statistical control in models (see Sect. 2.7). 

Since 1998, Observatoire PELAGIS has set up an internet site with contact informations 
and what to do when a stranded cetacean is found. Data on French household access to 
internet at the nation level were taken from "Centre de Recherche pour l'Etude et 
1' Observation des Conditions de Vie" 2 in order to incorporate this covariate as a statistical 
control in models (see 2.7 below). 

In this study, we aggregated the number of small delphinid strandings over a period of ten 
days at the county level. Thirty-three years of data from 1 980 till 201 2 were extracted for each 
county. Data are an array of three periods of 10 days times 12 months times 33 years times 
eight counties, that is a sample size of nearly 10,000 (3 x 12 x 33 x 8 = 9,504). Small 
delphinids are meant to include both the Common dolphins Delphis delphis and Striped 
dolphins Stenella coeruleoalba. These two species are similar in size, and were pooled in this 
analysis because carcasses were sometimes too putrefied to allow unambiguous species 
identification. Pooling makes the assumption that stranded carcasses of these two very similar 
species would trigger the same response by casual observers. 

Coastline lengths of sea-bordering counties were taken from the "Observatoire National 
de la Mer et du Littoral" . 3 Coastline lengths of the eight French Atlantic counties varies by 
one order of magnitude ("Pyrenees Atlantiques" < 100 km and "Finistere" > 1,000 km), 
and were incorporated in every model to account for the positive relationship between 
stranding probability and seashore length. 

Notations 

E(jc) and V(jc) are respectively the mean and variance of the random variable x. ^(X) 
denotes the Poisson distribution of mean parameter X. ^(".^(overdispersion, X) denotes the 
Negative Binomial distribution of mean parameter X and overdispersion parameter (Zheng 
et al. 2006). The usual parametrization of a Negative binomial (for example, Gelman et al. 

n 1 

2003 page 576) is Jf38{n, k) where X = — and overdispersion =1-1 — , but we chose to 

K K 



2 http://www.credoc.fr/. 

3 http://www.onml.fr/accueil/. 

Springer 



Biodivers Conserv 



follow the more convenient parametrization of Zheng et al. 2006. ^i(L,U) denotes the 
uniform distribution with lower bound L and upper bound U. J/~{p.-, cr) denotes the normal 
distribution of mean parameter /.( and scale parameter ct. ^(/i, a) denotes the Cauchy 
distribution of mean parameter fi and scale parameter a. ct, v) denotes the Student 
distribution of mean parameter /i, scale parameter a and v degrees of freedom, 
^(shape, rate) denotes the Gamma distribution parametrized in terms of a shape and rate 
(or inverse-scale) parameters. 



Statistical analysis of raw counts of stranded Cetaceans 



Strandings are count data, often analyzed with Generalized Linear Models or Generalized 
Additive Models with a Poisson likelihood (Pikesley et al. 2012; Fruet et al. 2012). The 
Poisson distribution is useful to model integer value variables (O ' Hara and Kotze 20 1 0), but it 
is restrictive because it assumes a variance equal to the mean (Berk and MacDonald 2007). 
Count data often exhibit two characteristics: (1) a variance larger than the mean, that is 
overdispersion (Berk and MacDonald 2007); and (2) an excess of zeros (Martin et al. 2005). 
The Negative Binomial distribution expands the Poisson distribution to account for over- 
dispersion (Berk and MacDonald 2007): the latter scales the variance with the mean 
(V(y) = overdispersion x E(y)). The Poisson distribution assumes overdispersion = 1 , 
while the Negative Binomial distribution assumes overdispersion > 1 . Overdispersion arises 
because observations are correlated (Berk and MacDonald 2007): in the case of marine 
mammals, the same process (for example epizootics (Carpenter 2013; Di Guardo 2012)) 
affecting a demographic unit may induce a multiple strandings event, thus violating the 
implicit assumption of independence between strandings under a Poisson likelihood. 

Excess zeros in count data may result from imperfect detection (Martin et al. 2005; Peron 
et al. 2010). Given the observational nature of strandings, detection probability may be an 
issue. Zero-inflated models offers a solution to model imperfect detection. For data y mod- 
elled with a Poisson likelihood, the frequency ratio of zero and single events is 
Pr(y = 0) 1 

— ; r = „, s - For data y modelled with a Negative Binomial likelihood, this ratio 

Pr(y = 1) Eft) 

Pr(y = 0) overdispersion ... . 

becomes — ; = — . Because the overdispersion parameter is at least 1, a 

Pr(y=l) E(y) P 

negative binomial distribution induces an excess of zeros compared to a Poisson distribution. 

The issue of choosing the appropriate likelihood is important to avoid overfitting (see 

Electronic Supplementary Materials). Overfitting is capitalizing on chance, and can happen 

when a statistical model is learning noise from a peculiar dataset rather than a meaningful 

signal that may be found in independent data (Babyak 2004). Overfitting can result in 

misleading inferences. 



Model building 



Let yijkt denotes the i datum corresponding the number of stranded carcasses to wash 
ashore during ten of days of month j and year t on the coastline of the k {h county. Month 
and year effects were always included in the variance components models we built. Month 
effects (at;) were assumed exchangeable and drawn from a normal distribution with 
common variance C7^ 0nth . Likewise year effects (/J,) were assumed exchangeable and 
drawn from a normal distribution with common variance (7? . To account for local 
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conditions at the county level, county effects (6^) were assumed exchangeable and drawn 
from a normal distribution with common variance cr^ ounty . 

Our simplest model accounting for the data structure is a cross-classified random 
intercept model (Browne et al. 2007), with month, year and county modelled as random 
effects. We were interested in assessing any temporal trend in the data. This trend may be 
slightly different for each county because of local factors not explicitly accounted for in the 
model. We thus included a county-level random slope, that is we modelled slope param- 
eters as exchangeable. This assumption is important: it is known a priori that "Les Lan- 
des" county differs from all the others owing to its institutionalized monitoring of beaches, 
and constant detection probability, since the 1990s. Violation of the exchangeability 
assumption of counties would betray non-constant detection probabilities in the other 
counties. The different parametric models we entertained are summarized in Table 1. 

Ml may be viewed as a default model: 

yj/*,~^(Coastjfc x l jkt ) (1) 

where the mean number (per 10 days and per unit of littoral length) of stranded carcasses in 
month j of year t found in county k is: 

\og{X jkt ) = 0 M + 9 k 2 x t + ctj + p, (2) 

M2 expands Ml to account for an excess of zeros. The likelihood for the data is a mixture of a 
degenerate mass at zero and a Poisson distribution. The model for the mixing proportion is 
also a cross-classified random intercept model as described above. M2 has approximately 
twice as many parameters asMl.M3 is a generalization of Ml and includes an overdispersion 
parameter common across all years. M4 further generalizes M3: the overdispersion param- 
eters are year-specific. M5 is an elaboration of M4 described below (see Sect. 2.6). Models 
Ml— 2 assumed a Poisson likelihood for the data, while models M3— 5 assumed a Negative 
Binomial distribution. In all models, a logarithmic link function was used. In zero-inflated 
models, the mixing proportion was modelled with a logit link. 



Structuring overdispersion 

Overdispersion is usually treated as a nuisance parameter, or a parameter not of primary 
interest itself, but that must be taken into account for correct inferences. Yet, overdis- 
persion may capture interesting information (Zheng et al. 2006). Our large dataset enables 
us to estimate a year-specific overdispersion parameter (Model M4, see Tables 1 and 2). In 
M4, overdispersion parameters are unstructured: we can obtain a time-series of overdis- 
persion parameters over the study period. The statistical problem becomes to identify 
'surprising' points. Fuquene et al. (2014) proposed a method to identify outliers and 

structural breaks in time-series. Let k, = , so that k, > 0. The k, > i are 

overdispersion, — 1 

modelled with a random walk, on a logarithmic scale to ensure positive values: 

{log(Kt) + e t 
log(Ki) + J>i 
i=l 

The next value of a given log(jc t +i) is the previous value log(jc f ) plus a small deviation e t . 
Large deviations in magnitude may betray outliers. Fuquene et al. (2014) proposed to 
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assume a heavy-tailed distribution, the Student t, for the £,. The Student t distribution can 
be represented as a scale mixture of normal (Gaussian) distributions (Andrews and Mal- 
lows 1974): let z, be a random variable with a normal distribution, z, ~ Jt (0, cr £ ); and co, 

an independent, positive, random variable. The random variable £, = Z ' is a scale 

mixture of normals. The conditional distribution of e, given co, is ^"(0,— — ). If 

co f ~^(|,|), the marginal distribution of e, is Student r with v degree of freedom. This 
representation of the Student t as a mixture of normals with different variances is useful: 
the parameter co, has mean 1, the conditional distribution of e, is normal with scale ct £ . 

When co, < 1 , > 1 allows for a greater variance than usual: the deviation e, may have 

a larger than usual magnitude. Thus a value of co, < 1 will betray an outlier in the time- 
series. The parameter v for degrees of freedom controls tail-heaviness (or how far the tail 
extends from the mean value): for large values of v (v> 30), the Student t distribution is 
similar to a normal distribution, while small values of v authorized outliers. We fixed v = 4 
as in Fuquene et al. (2014). Fuquene et al. (2014) further proposed a weakly-informative 
Beta2 prior for the scale parameter ct 6 . We embedded Eq. 3 in model MA (Table 1) to 
detect outliers in the time-series of overdispersion coefficients. This model is labeled M5 
(Table 1). 

Confounding variables 

A confounder is a variable known to be causally related to the response variable, and 
whose omission in an analysis would generate a spurious relationship between the latter 
and a third variable of interest. We distinguished between confounding and lurking vari- 
ables which are respectively known and unknown to the modeller. Known confounders in 
our case are (1) demographic evolution at the county level, (2) household internet access 
and (3) law enforcement. The more populous a county, the more likely a stranded carcasses 
may be reported. A 1988 law enforced mandatory reporting of cetacean strandings to 
Observatoire PELAGIS and RNE by town and village administrations. The passing of this 
law may have triggered an increase in reporting rates from local authorities. In 1998, 
Observatoire PELAGIS set up an internet site with information on whom to contact and 
what to do after the discovery of a stranded animal. Reporting rates from citizens may 
consequently have increased with household internet access. Upon selecting a model with 
an appropriate likelihood, we included these variables in Eq. 2: 

\og(lj kt ) = 8 k j + Q ki 2 x t + 8 k3 x law, + 6 kA x internet, + 9 kt5 x pop size t , + a,- + f}, 

(4) 

where law, = 0 if t< 1989 and 1 otherwise; internet, denotes the percent of household with 
internet access from 1998 onwards; and pop size^. , is the population size of county k in 
year t (see Sect. 2.2). Positive effects of all three covariates are expected if they are true 
confounders, in which case we also expect the regression coefficient for time t to attenuate 
or zero-out. Finally, we also compared the fit (via the EarthMover distance, see Sect. 2.8) 
of a model in which the parameters 9 k ^, that is the regression coefficients for time, are set 
to zero in Eq. 4 to the fit of our parametric model MA (Table 1). 
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Model checking, selection and testing 

For each model, we used 3,000 draws from the posterior to predict each datum, and 
computed the average histogram of predictions ("posterior predictive checks", Gelman 
et al. 1996). This histogram was compared to the histogram of observed data, and model fit 
was visually assessed (see Fig. 3 in Minami et al. 2007 for a similar exercise). To 
numerically quantify model fit, we used the EarthMover distance which has an intuitive 
interpretation for the comparison of two histograms: each histogram may be viewed a pile 
of dirt and the EarthMover distance reflects the amount of dirt times the distance needed to 
turn one pile into the other. Each datum is like a sand grain, and the aim is to quantify how 
many sand grains and how far do we need to move them to transform the histogram of 
predictions into the empirical histogram of observed values. For inference, we selected the 
model with the smallest EarthMover distance, that is the model which was able to generate 
a histogram of predictions the most similar to observations. 

Model fitting and priors 

We favored a Bayesian framework. Bayesian estimation of random effects with a small 
number of levels (12 months, eight counties) is adequate (Stegmueller 2013). We used 
STAN (Stan Development Team 2013) called from R (R Development Core Team 2013) 
with the package rstan (http://mc-stan.org/). Weakly informative priors were used 
(Gelman 2006; Gelman et al. 2008). For variance parameters, we used a half-Cauchy prior: 
a~-^ + (0, 1). For intercept parameters, we used a Student t: intercept ~ .5^(0, 10,7). For 
slope parameters, we used a Student t: slope ~ .9^(0, 2.5, 7). For overdispersion parameters 
in models M3 — 4, we used uniform priors: overdispersion ~ 100). Three chains were 
initialized with different starting values. After appropriate burn-in (1,000 iterations) and 
thinning of the chains (1 value every two iterations stored), convergence was assessed 
using the Gelman-Rubin convergence diagnostic (Cowles and Carlin 1996). Unless stated 
otherwise, posterior means are reported with 95 % Highest Probability Density intervals 
( 2 .5%Mean 9 7 5%) (Louis and Zeger 2009). Inferences are based on a posterior sample of 
3,000 iterations. Covariates were standardized following (Gelman 2008). STAN and R 
codes are available as supplementary materials. 

Results 

The mean number of stranded small-sized delphinids per 10 days over the study period was 
lEfjobs) = 1 an d the variance was V(y 0 b s ) = 34. Seventy-three percent (73 %) of the data 
were zeros. All models converged (r < 1.1). Models with a Poisson likelihood did not fit 
the data well: these models failed to predict enough zeros and predicted too many 
observations between 1 and 10 strandings (Fig. 2a). Changing the likelihood for a Zero- 
Inflated Poisson fixed the problem for the excess zeros, at the cost of twice as many 
parameters; but the models still overpredicted frequencies of strandings of 1-10 individuals 
(Fig. 2b). Using a Negative Binomial likelihood however resulted in a more acceptable fit 
(Fig. 2c, d). This ranking is similar to that of Minami et al. (2007) with shark bycatch data. 
Poisson and Zero-Inflated Poisson models were unable to predict events of more than 50 
stranded dead dolphins. Negative binomial models were able to predict events of up to 80 
stranded dead dolphins. No models could satisfactorily take into account the handfull of 
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Fig. 2 Parametric model fit. After model fitting, the posterior distribution of estimated parameters was used 
to generate replications of the data. For each replicate, the empirical histogram of predictions was computed. 
The mean predicted histogram were finally compared to the target histogram of observed data to assess 
whether estimated parameters were able to predict similar data sets to the one the model was calibrated with. 
Both x- and y-axes are on a square-root scale for ease-of-read purposes. Histograms are truncated at 100 
(only five events, among 9,504, of more than 100 stranded carcasses over 10 days during the study period) 

extreme events: five events of > 100 stranded carcasses, including one with 349 animals. 
Model fit is numerically summarized in Table 2. The best model was MA: a Negative 
Binomial likelihood provided an aquedate fit to these data. The large sample size of our 
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Fig. 3 Model fit at the level of individual counties. We compared within each county predicted and 
observed histograms. The mean prediction is depicted in solid line, and a 95 % HPD interval as dashed lines. 
Both x- and y-axes are on a square-root scale for ease-of-read purposes. Histograms are truncated at 100 
(only five events, among 9,504, of more than 100 stranded carcasses over 10 days during the study period) 



data allowed us to estimate a year specific overdispersion parameter, which also slightly 
improved model fit (Table 2). Model fit at the county level was also adequate (Fig. 3). 

Figure 4a shows the ranking of each factor (month, year or county) in relative impor- 
tance: monthly variation were the most pronounced in magnitude, followed by county- 
level variations, and lastly by year-level variations. These results are in line with those of 
Peltier et al. (2014) who used a different modelling approach with a subset of the data we 
used. Estimated parameters can be interpreted as multiplicative factors (Fig. 4b, c): the 
number of reported stranded carcasses of small delphinids increased 5-fold in February 
compared to an average month like May or December, or a 10-fold increase compared to a 
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Fig. 4 Visualizing results, a Ranking of sources of variation in terms of importance and magnitude. 
Seasonal and yearly variations were most and least pronounced respectively, b Estimated multiplicative 
factors for monthly variations, c Estimated multiplicative factors for yearly variations. The year of law 
enforcement for mandatory stranding reporting is indicated with an arrow. Estimation uncertainty is 
displayed with shading using the R package denstrip (Jackson 2008) 

summer month like July or August (Fig. 4b). Strandings were less frequent during the 
summer and autumn month but increased dramatically during winter and spring months 
(Fig. 4b, Peltier et al. 2012). Yearly variations were present but less pronounced, except for 
a spike in the late 1980s. This spike happened the same year a law for mandatory reporting 
of stranded animals was enforced. 

Between 1980 and 2012, the mean number of strandings increased 2.37.O12.3 fold over 
the Bay of Biscay coastal areas (Fig. 5i). This overall trend did not mask county-level 
peculiarities: small delphinid strandings increased in all the 8 French counties bordering 
the Atlantic Ocean (Fig. 5a-h). Using Ml, a model that did not fit the data (see Table 2), 
would have resulted in different inferences regarding trends (see Supplementary Fig. SI). 
The model we selected included overdispersion parameters (Fig. 6). Overdispersion, that is 
a variance inflating factor relative to a Poisson distribution, decreased over the 33 studied 
years (Fig. 6a). We further applied the method described in Fuquene et al. (2014) to 
identify structural breaks in this time-series (Fig. 6b, c). The years 1992 and 1998 were 
outliers, both because overdispersion decreased sharply compared to the rest of the series. 

Discussion 

We analyzed our data at the county level to capitalize on known differences between 
counties. In particular, monitoring has been institutionalized since 1989 in "Les Landes" 
county. Local permanent staff is in charge of monitoring and cleaning the beaches that 
borders the Bay of Biscay for aesthetic reasons linked to tourism. We can reasonably 
expect reporting rates of stranded cetaceans to have remained fairly constant since the 
1990s there. Yet, the overall shape of the trend in "Les Landes" county does not differ 
from the other counties (Figs. 5, 7). "Les Landes" county has an overall higher number of 
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Fig. 5 Temporal trends in small delphinid strandings. Estimated trends from model M4 are depicted for 
each county (a-h), and the Bay of Biscay as a whole (i). The y-axis (on a logarithmic scale) represents the 
estimated mean number (per 10 days and per 100 km of coastline) of stranded carcasses, that is log(Aju) = 
Sit i + 8k2 x r is plotted for each county k. Dots are the raw data, solid line is the posterior mean, and dashed 
lines represent a 95 % HPD interval around the mean (without taking into account the year random effects 
fl n see Eq. 2 and Fig. 4c). In each of the county bordering the Bay of Biscay, the trend is that of an increased 
number of stranded small delphinids. Although counties have different intercepts, slopes over the study 
period are similar, suggesting no major difference in the processes behind these trends between the different 
counties. This pattern suggests that reporting rates, which are known to have been constant in "Les Landes" 
(g) since the 1990s, have also remained constant in the seven other counties 
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Fig. 6 Detecting outliers, a 
Time series of estimated 
overdispersion parameters from 
semi-parametric model M4. the 
inverse of overdispersion is 
plotted to constrain values in the 
unit interval [0:1], and avoid 
visual distortions because 
overdispersion coefficients are 
bounded between 1 and +oo. 
Dots represent the mean, and 
arrows a 95 % HPD interval, b 
Time-series of year-specific 
deviations (see Sect. 2.6). The 
solid line represents the mean, 
and dashed lines a 95 % HPD 
interval, c Identifying outliers 
and structural breaks when the 
parameter to, < 1 (see Sect. 2.6. 
Dots represent the mean, and 
arrows a 80 % HPD interval 



strandings (Fig. 7). This results from a favorable substrate and geography, 105 km of 
gentle-sloped sand beaches, which makes it a popular domestic touristic destination; and 
the absence of a major estuary. "Gironde", "Charente Maritime", and "Loire-Atlantique" 
counties take their name from large rivers that flow into the Bay of Biscay. "Pyrenees 
Atlantiques", "Morbihan" and "Finistere" counties have rocky and cliffy seashores, which 
lower the probability that a cetacean will strand and be detected. However, the estimated 
slope for "Les Landes" counties was comparable to that of two other counties ("Vendee" 
and "Pyrenees Atlantiques"), which vindicated the exchangeability assumption. By ana- 
lyzing our data at the county level, we were interested in comparing "Les Landes" to the 
other counties, and in particular to the "Gironde" county with which it shares many 
features (same substrate, tourism, geographical proximity). These two counties have 
similar trends over the study period, which further suggests that detection probability is not 
sufficient to explain away the observed increase in stranding numbers. 

Potential confounders 

The number of stranded carcasses is the result of many underlying processes (Peltier et al. 
2012): a change in any component may trigger a trend while other components remain 
unchanged. The observational nature of stranding data limits the scope for causal infer- 
ence. In particular, detection probability of stranded carcasses is an obvious confounder. 
To account for this, we incorporated several confounders in our final model (Table 3). In 
the final model with no covariate but 'time', the slope coefficient was estimated at 
o.66l .442.23- This coefficient decreased to -0.23O.962.24 when confounders were included, 
suggesting that the increasing trend may partly result from an increase in detection. 
However, the sign of the coefficients for household internet access and county population 
size have the wrong sign (Table 3). For the latter, the negative sign implied that more 
populous county had fewer reported strandings, which does not make sense if detection 
probability is related to the number of inhabitants; unless this is also correlated with other 
cultural changes by which the lay public would pay less attention to natural events. The 
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Fig. 7 Estimated county-level intercept (flji) and slope (B k2 ) from model M4. The .v-axis represents the 
county-specific intercept, and the y-axis the county-specific slope (for 'time'). "Les Landes" county has a 
larger intercept, that is a higher number of strandings than any other other county. Its slope coefficient for 
'time' is however comparable to other counties. The black dot symbolizes the grand mean for the Bay of 
Biscay 



Table 1 Parametric models: five models of varying complexity were considered; see Sect. 2.5 



Name Likelihood 



Mean function 



Overdispersion 



A/1 

M2 



Mi 
MA 
M5 



Poisson 

Zero-Inflated 
Poisson 



Negative Binomial 
Negative Binomial 
Negative Binomial 



log(% fa ) = () u + 0 k2 x r + Oj + P, 
>-jk, = njkr x 0+ 

(1 - n jkt ) x exp(O kA + 0 k2 x t + a. jA + /i, ,) 
where n jk , = logir'^o + i j2 + ji tl ) 
log(V) = 6 kA + 0 k2 x t + OLj + p, 
log(X jkt ) = 0 kA + 0 k2 x ( + aj + p, 

l0g(V) = °K\ + 0 k2 X 1 + CLj + P, 



overdispersion 
overdispersion, 

overdispersion, — 1 H 

K, 

K\ > 0 
for r> 1, 
log((C, + i) = log(K f ) + E, 



negative sign of the regression coefficient for household internet access is likewise sur- 
prising. In the era of omnipresent new technologies, it is paradoxical that a greater access 
to internet access should result in less reporting assuming that outdoor recreation habits of 
citizens have not changed. How can we explained the negative regression coefficient then? 
Household internet access is strongly correlated with 'time' over the study period, as can 
be seen from the large negative correlation between the regression coefficients for these 
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Table 2 Model selection with the EarthMover distance: the model which was able to make predictions 
most similar to observed data was selected 



Model EarthMover distance 



Ml 0.443 

M2 0.384 

M3 0.203 

M4 0.198 



The most richly parametrized model, M2, was not the one that provided the best fit because the likelihood 
was inappropriate for the small delphinid stranding data in the Bay of Biscay. Models with a Negative 
Binomial likelihood were clearly superior 



Table 3 Estimating the effect of potential confounders: estimated regression coefficients of potential 
confounders affecting detection probability are reported on the diagonal 



Covariate Time Law internet Pop size 



time 


-0.23 0.96 2.24 -0.80 


-0.86 


-0.10 


law 


-0.14 0.64 1.36 


0.63 


0.01 


internet 




-2.04 -0.40 1.10 


0.03 


pop size 






-2.37 -0.92 0.37 



Off diagonal figures correspond to empirical posterior correlations between parameters, and indicate col- 
linearity between 'time' and household internet access 



two covariates (Table 3): any increase in the regression coefficient for 'time' goes hand in 
hand with a similar decrease in the regression coefficient for household internet access. 
One solution to resolve this collinearity problem is to delete one of the covariates. When 
the regression coefficient for time was set to 0 in Eq. 4, model fit was slightly worse 
(EarthMover distance = 0.208 versus 0.198 for a model with just 'time'). In this case, the 
sign of the regression coefficient for household internet access became positive 
(o.060.75i.57), that of law enforcement remained positive (0.69LI31.59) and that of popula- 
tion size remained negative (-2.32 — 0.84o.44). 

Law enforcement in 1988 probably had a causal impact and increased detection 
probability: we observed a positive regression coefficient for this indicator variable. We 
detected a structural break in the time-series of overdispersion parameters corresponding to 
the year 1992 (Fig. 6c): overdispersion decreased in that year. Passing the law for man- 
datory reporting of stranded carcasses by local French administrations triggered an increase 
in the mean number of reported stranding (Fig. 4c), but compliance may have been partial 
for the first couples of years after law enforcement. A period of partial compliance could 
reflect the time needed for local authorities to adapt and implement routine procedures 
when a stranded carcass is discovered. By 1992 however, full compliance was probably 
achieved, which may account for the observed decrease in variability. Another, non- 
mutually exclusive factor, may be a greater public awareness of dolphin strandings fol- 
lowing the Mediterranean morbilivirus epizootics of the early 1990s (Aguilar and Raga 
1993), which was prominently covered by both local and national media. 

Evidence for an effect of household internet access on reporting a stranded animal to 
Observatoire PELAGIS is less clear. During the early 1990s, less than 5 % of French 
household had an internet access and Observatoire PELAGIS had no website with instruc- 
tions to follow in case of discovery of a stranded carcass. Yet, the overall and local trends were 
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all one of increase during the 1990s (Fig. 5). If internet access had a causal impact of reporting 
rate, it would manifest itself from the 2000s onwards, but this period corresponds to reduced 
rate of increase compared to the preceding decade. The internet covariate is 0 or very small for 
most of the 1990s. It thus appears unlikely that access to new technologies had an effect on 
detection in the 1990s which was nevertheless a period of increase in strandings. This can 
explain why the fit of model without 'time' is slighly worse than one with 'time'. Finally, 
although population size on the French Atlantic coast increased over the study period, the sign 
of the estimated coefficient for this confounding variable is not coherent with a causal effect 
on detection probability. This in turn suggests a preponderant role of local administrations 
and their staff (e.g. earthmovers) in reporting stranded carcasses, over that of the lay public. 

Public awareness of cetacean conservation issues has likely increased over the study 
period. Were the lay public the only channel through which strandings were reported to 
Observatoire PELAGIS, population size would have increased reporting rates in all 
counties but "Les Landes" . The reason why this does not show up in the data is probably 
linked to the effectiveness of the 1988 law. Most reportings are done through the "official 
channel" of municipalities calling directly Observatoire PELAGIS. Although Observatoire 
PELAGIS does get more calls from the lay public, these calls are largely redundant with 
those of municipal staff. This might be the reason why reporting rates may have remained 
more or less constant since the 1990s. Population size may also be too crude a proxy: more 
relevant information could be changes in leisure activities, either rural or urban. The 
former is more in favor of higher reporting rates than the latter, yet the French population is 
currently more urban than ever before. 

Causal inference with observational data is inherently more difficult and tentative than 
with experimental data (Goldthrope 2001): without experimental control, an unknown 
confounder may still be lurking. Here, we focused on whether observed trends could have 
arisen from variations in reporting rates alone. We concluded that this may have been the 
case in the 1980s, but unlikely afterwards. 

What happened in the 1990s? 

A major threat to marine mammal populations worldwide is bycatch by fisheries targeting 
commercial fish species (Morizur et al. 1999; DeMaster et al. 2001; ICES 2005; Read et al. 
2006; Moore and Read 2008; Reeves et al. 2013; Thompson et al. 2013). The Bay of 
Biscay corresponds to ICES Ecoregion G, subdivisions Villa and VHIb. In the 1980s, the 
number of pair trawlers operating in the whole ICES Ecoregion G Subdivision VIII 
increased from less than 50 boats to more than 150 in the early 1990s and down to 100 in 
the early 2000s (Villalobos Hortiz 2008). The number of purse seiners remained more or 
less the same (between w25 and «40) during the 1990s but doubled in the early 2000s 
(Villalobos Hortiz 2008). Dolphin bycatch by pair trawlers or purse seiners is documented 
(Read et al. 2006; de Boer et al. 2012; Thompson et al. 2013), and a represent a serious 
threat to common dolphin population viability on the Bay of Biscay (Mannocci et al. 
2012). In particular, common dolphins are at risk with the sea bass trawl fishery operating 
the Bay of Biscay because of the dietary overlap between sea bass (Dicentrarchus labrax) 
and common dolphins (Spitz et al. 2013). Because of high energy requirements (Spitz et al. 
2012), small delphinids target preferentially prey of high energy contents such as European 
anchovies (Engraulis encrasicolus) and sardines (Sardina pilchardus) (Meynier et al. 
2008; Spitz and Jouma'a 2013), on which European sea bass also feeds. A byproduct of 
dietary overlap may be increased interactions between dolphins and commercial pair 
trawlers, and eventually dolphin bycatch (Mannocci et al. 2012; Spitz et al. 2013). 
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We acknowledge that inferring the causes behind the trend observed in our data is a 
difficult endeavour, and ruling out potential confounders (as in Jepson et al. 2013 or Wright 
et al. 20 1 3) is certainly a more realistic goal than finding the unique causal mechanism behind 
the observed pattern. Yet, there is ample independent data on the threats that bycatch poses to 
marine mammals (Morizur et al. 1999; DeMaster et al. 2001; ICES 2005; Read et al. 2006; 
Moore and Read 2008 ; Reeves et al. 20 1 3 ; Thompson et al. 20 1 3) to argue for at least a partial 
role of the later in the average 7-fold increase in stranding numbers of small delphinids we 
observed in the Bay of Biscay over 30 years (Fig. 5). Excluding data from the 1980s, the 
average increase over the whole Bay of Biscay is 2.03.85.7-fold. Using an explicit physical 
drift model to account for stranding probabilities of dead dolphins, Peltier et al. (2014) 
demonstrated that the continental shelf of the Bay of Biscay is a mortality hotspot for common 
dolphins. Peltier et al. (2014) further estimated that a 2-fold increase in stranded common 
dolphins between winter and summer was due to drifting conditions alone. Here, we esti- 
mated a 10-fold increase (Fig. 4b), that is a 5-fold increased compared to what would be 
expected with seasonal drift conditions alone, which is again similar to what Peltier et al. 
(2014) found using a different approach that accounts explicitly for drift. 

We discussed earlier a structural break we identified in 1992, but set aside the other one 
identified in 1998. Close inspection of Fig. 6 and data suggest that the outlier is rather the 
year 1997, in which overdispersion increased sharply compared to the surrounding years. 
This increase was due to the occurrence of six extreme events of more than 50 strandings 
per 10 days in mid-February 1997 in five counties. Examination of individual cases by 
querying our database revealed that a majority of the animals showed lesions either 
characteristic or suggestive of interactions with fishing gear (Kuiken 1994). Although our 
selected model did not cope well with predicting extreme stranding events, in this case the 
mean trend was unaffected (Fig. 5), and the overdispersion parameter acted as a buffer 
against these outlying observations. This, in turn, shed light on why a Negative Binomial 
likelihood is adequate for small delphinids data: multiple stranding events can stem from 
bycatch of a large number of animals that were interacting with a commercial fishery. The 
stranding events are no longer independent, but are correlated through a common cause of 
accidental death. Our study thus suggests using overdispersion as an indicator of years of 
unusual mortality. Yet, the precise identification of underlying causes behind unusual 
events is not something our models can elucidate with certainty, and ancillary information 
will be needed. In other words, overdispersion as an indicator has a low specificity, but this 
is unsurprising given the phenomenological nature of the model. In fact, the year 1997 is 
also an outlier in that several winter storms resulted in favorable drift conditions, which 
may partly explained the surge in observed strandings. However, Peltier et al. (2014) also 
found that 1997 was the first year with anomalously high strandings even after accounting 
for drift conditions. In our case, in spite of drift being a potential confounder, it does not 
alter how a large proportion of stranded small delphinids in the beginning of 1997 showed 
lesions characteristic of bycatch. The merit of model M5 is to allow robust (with respect to 
noise) inference on the mean trend, to provide a parameter that can be linked to a common 
cause of mortality (bycatch, epizootics), and to suggest in a principled way what are 
surprising observations in the time-series of these parameters. 

Conclusion 

We have provided an in-depth analysis and description of small delphinid strandings in the 
Bay of Biscay, capitalizing on county-level specificities to critically assess the consistency 
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and quality of data collected by the French stranding network RNE since 1980. With no 
claims of novelty (for example, see Minami et al. 2007), we suggested to use variance 
component models to take into account some peculiarities of count data, and proposed 
biological interpretation of parameters of interest. A good fit with the observed data was 
obtained with relatively simple models including seasonal and yearly variations and a trend 
over time. The model allowed to rank sources of variations in importance (Fig. 4a), and to 
reduce a large data set of nearly 10,000 points to less than 100 interpretable parameters. 
The number of stranded small delphinids have increased seven-fold since the 1980s. We 
discussed causal inference with observational data and investigated potential confounders 
related to detection of stranded carcasses. We inferred that the French stranding records 
appear most homogeneous in effort since the 1990s, with the corollary that the observed 
increased in standings in due to other factors than an increase in reporting rates. 

Although we adopted a narrow focus on the Bay of Biscay and RNE in order to provide 
a documented, context specific, analysis of observational data, we feel that our models may 
be of interest to other national stranding networks for monitoring purposes. Moreover, the 
graphical representation of results is relatively new and illuminating for applied research 
(Gelman 2005; Hector et al. 201 1). Recent research also suggests that traditional reporting 
of results by means of F-statistics may create an illusion of predictability, even for experts 
(Soyer and Hogarth 2012). 

We discussed the choice of an appropriate likelihood for the data, which seems an 
overlooked issue in the marine mammal stranding literature, and found the Negative 
Binomial useful and interpretable in the context of small delphinid standings. This need 
not be the case with other taxa, especially species which are less abundant (baleen whales 
for example), but we suggest it should be tested along with the more common Poisson 
model to avoid overfitting (see Online Appendix). Finally, we did not consider modelling 
drift explicitly as in Peltier et al. 2013, 2014. The latter is to be preferred since any signal 
can be more meaningfully interpreted as freed from the confounding influence of drift. 
However, our aims here were to analyze standings data in a general setting (see Peltier and 
Ridoux 2013 for setting the scene of stranding data analyses with a drift model). A future 
extension is to merge this work with explicit drift modelling (Peltier et al. 2013, 2014) for 
assessing and monitoring Good Environmental Status with respect to marine mammals 
within the EU Marine Strategy Framework Directive. 
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